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The density and temperature dependence of the nuclear symmetry free energy is investigated using 
microscopic two- and three-body nuclear potentials constructed from chiral effective field theory. The 
nuclear force models and many-body methods are benchmarked to properties of isospin-symmetric 
nuclear matter in the vicinity of the saturation density as well as the virial expansion of the neutron 
matter equation of state at low fugacities. The free energy per particle of isospin-asymmetric nuclear 
matter is calculated assuming a quadratic dependence of the interaction contributions on the isospin 
asymmetry. The spinodal instability at subnuclear densities is examined in detail. 

I. INTRODUCTION 

Determining the thermodynamic equation of state (EoS) of nuclear matter is a central objective in 
modern nuclear theory. The isospin-asymmetry dependence of the EoS is essential for many phenomena 
in nuclear physics and astrophysics [T]. Next-generation radioactive beam facilities Hi studying the 
reactions and structure of exotic neutron-rich isotopes in particular provide motivation to improve our 
microscopic description of highly isospin-asymmetric nuclear matter. To a certain extent, the isospin- 
asymmetry dependence of the EoS is described by the so-called symmetry free energy. In this work, 
starting from microscopic calculations of the EoS of isospin-symmetric nuclear matter and pure neutron 
matter using two- and three-body chiral nuclear interactions, we examine in detail the density and 
temperature dependence of the symmetry free energy. Furthermore, we construct the EoS of isospin- 
asymmetric nuclear matter, and study the behavior of the nuclear liquid-gas instability as the proton 
fraction is decreased. The present work is a first step toward the development of a chiral effective 
field theory thermodynamic equation of state across the temperatures, densities and isospin asymmetries 
relevant for describing astrophysical phenomena and the matter produced experimentally in heavy-ion 
collisions at moderate energies. 

Chiral effective field theory (yEFT) provides the basis for the study of strongly-interacting matter 
at the energy scales characteristic of normal nuclei gHS]. In yEFT, microscopic nuclear interactions 
are organized in a systematic expansion, with many-nucleon forces naturally included. The low-energy 
constants parametrizing the interactions are generally fixed by high-precision fits to nucleon-nucleon 
scattering phase shifts and properties of light nuclei. Employing chiral interactions in calculations of 
nuclear many-body systems then gives pure predictions without additional fine tuning, and theoretical 
uncertainties can be estimated [BlUO] by varying the resolution scale, the fitting procedures applied 
to fix the low-energy constants, and the chiral order of the nuclear potentials. In the case of isospin- 
symmetric nuclear matter, semiempirical constraints from the zero-temperature saturation energy, density 
and incompressibility as well as the critical point of the nuclear liquid-gas phase transition have been 
reproduced with low-momentum chiral nuclear forces in many-body perturbation theory [111 112] . This 
motivates a study of isospin-asymmetric nuclear matter, which is by comparison much less constrained 
by experimental data. 

The symmetry free energy Fsyra{T,p) is defined as the difference between the free energy per particle 
in homogeneous isospin-symmetric matter (SNM) and pure neutron matter (PNM): 

Esym(r, p) = F[T, P,5=l)- F{T, P, 5 = 0), (1) 

where T is the temperature, p = pn + Pp is the total nucleon density, and 5 = (pn — Pp)/(Pn + Pp) is the 
isospin-asymmetry parameter (with p^/p the neutron/proton density). The free energy per particle of 
homogeneous nuclear matter with proton fraction Tp = (1 — 5)/2 can be written as 

F(T, p, J) = F{T, p, (5 = 0) + F,y„,(T, p) (3(T, p, J) <5^, (2) 

where (3(T, p, J = 1) = 1. If isospin-symmetry breaking effects are neglected (3(T, p, 5) is an even function 
of 5. It has been validated in various microscopic many-body calculations (see e.g.. Refs. mm) that 
the isospin-asymmetry dependence of F{T,p,6) at zero temperature is approximately quadratic to high 
accuracy over the entire range 0 < J < 1, i.e., (3(T = 0, p, S) ~ 1. At finite temperatures however the free 
Fermi gas contribution to F{T,p,S) contains large terms with quartic and higher powers of 6, which we 
quantify explicitly in this work. By comparison, in initial calculations we have found that at the densities 
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and temperatures relevant for the nuclear liquid-gas phase transition the interaction contributions give 
rise to weaker nonquadratic terms and will therefore be assumed to have a quadratic dependence on 6 in 
this work. Future research will address the accuracy of this approximation in greater detail. 

The liquid-gas phase transition in isospin-asymmetric nuclear matter (ANM) involves isospin distilla¬ 
tion: in the transition region the system separates into two phases whose proton concentrations deviate 
from the global value Fp, with 0 < < Yp and Yp < < 0.5 for the case Yp < 0.5. These 

distillation effects are a generic property of first-order phase transitions in binary thermodynamic sys¬ 
tems. If isospin-symmetry breaking effects are neglected, neutrons and protons are thermodynamically 
indistinguishable in SNM. Hence, thermodynamically SNM is a pure substance, and there is no isospin 
distillation for <5 = 0. The new features of the phase transition in ANM were discussed in detail in Refs. 

using different phenomenological models of the nuclear force. The transition region is comprised of 
regions of metastable and unstable single-phase equilibrium, corresponding to different dynamical phase 
separation mechanisms: nucleation and spinodal decomposition [ 101121 ]. In this work we focus mostly on 
the spinodal which delineates the inner region of thermodynamic instability where no metastable state 
can exist. The evolution of the unstable spinodal region with increasing isospin asymmetry is analyzed 
in terms of the trajectory of the critical temperature Tc{d). Moreover, we determine the neutron drip 
point in cold nuclear matter and the fragmentation temperature above which no self-bound drop of liquid 
nuclear matter can exist. It should be emphasized that in the present paper we discuss the liquid-gas 
instability of infinite nuclear matter, i.e., bulk nucleonic matter without Coulomb interactions. The in¬ 
clusion of surface energies and the Coulomb repulsion of protons is required for an accurate description 
of the matter produced in intermediate-energy heavy-ion collisions. In neutron stars and core-collapse 
supernovae the realization of the nuclear liquid-gas instability is strongly affected by the presence of 
a (highly incompressible) charge neutralizing background of electrons (and myons) [221 123]. In partic¬ 
ular, the competition between nuclear and Coulomb interactions (frustration) entails the fomration of 
mesoscopic inhomogeneities with nontrivial spatial structures. These so-called pasta phases have been 
studied extensively in the literature [IMTI] . Including these effects as well as the presence of few-nucleon 
bound-states [32H84] at very low densities represents a future challenge. 

The paper is organized as follows. In Sec. [IT] we recall the main results for the EoS of SNM obtained in 
Ref. [12] . and show results for additional derived thermodynamic quantities, i.e., the entropy per nucleon 
and the internal energy per nucleon. In Sec. |III| we extend the calculations to pure neutron matter. The 
zero-temperature results are compared to those from recent quantum Monte Carlo simulations while the 
finite-temperature EoS at low densities is compared to the virial expansion. In Sec. |IV| we investigate 
the temperature and density dependence of the symmetry free energy, entropy, and internal energy. The 
thermodynamics of isospin-asymmetric nuclear matter is studied in Sec. [^ In particular, we examine in 
detail the dependence on isospin asymmetry of the EoS of a free nucleon gas. Finally, Sec. jVl] provides a 
short summary. 


II. ISOSPIN-SYMMETRIC NUCLEAR MATTER 

In Ref. [H] we calculated the free energy per nucleon in infinite homogeneous SNM using the Kohn- 
Luttinger-Ward |3S1[5^ many-body perturbation series including contributions up to second order, i.e., 

F{T, p,S = 0)= FoiT, fio) + F,ei{T, po) + AFi(T, po) + po) + 0{X^). (3) 


Here, A counts the number of interaction insertions, Fq(T, pq) corresponds to a nonrelativistic free nucleon 
gas, and Fi.ei(T, /tq) is a correction term which together with Fq reproduces the properties of a relativistic 
free nucleon gas over a wide range of densities and temperatures m- The first- and second-order terms 
Fi{T,plq) and F 2 {T,pq) receive contributions from both the two-body and the three-body nuclear force. 
The second-order term F 2 {T, po) includes (temperature and density dependent) self-energy corrections, 
and has been evaluated by approximating the three-nucleon interaction with a temperature and density 
dependent effective two-body potential (for details see Refs. [T2l[38l - I^ . Explicit formulas for the different 
contributions in Eq. (§ are given in Ref. m- The effective one-body chemical potential /ig is in one-to- 
one correspondence with the nucleon density via 


P{T, Po) 


/dk 


1 -b exp 


k^/2M - po 


( 4 ) 


where r G { — 1/2,1/2} is the isospin projection quantum number and M ~ 938.9 MeV is the average 
nucleon mass. For Eq. (§ to be sufficiently converged at second order in A, low-momentum interactions 





3 


have to be used, i.e., interactions with restricted resolution in coordinate space (corresponding to an ul¬ 
traviolet cutoff in momentum space). The various sets of N3LO (i.e., fourth order in the chiral expansion) 
two-body and N2LO three-body chiral low-momentum interactions used in Ref. m correspond to differ¬ 
ent regularization methods, resolution scales A, and low-energy constants. For interactions constructed 
at resolution scales A < 450 MeV appropriate perturbative behavior was found. The SNM equation of 
state obtained from the sets of two- and three-body potentials denoted by n31o414 {A = 414 MeV) and 
n31o450 {A = 450 MeV), respectively, (see Refs, [m WA HH] for details) agree with empirical constraints 
from the zero-temperature saturation energy, density and incompressibility [44H47j . and with estimates 
for the critical point of the nuclear liquid-gas phase transition obtained through the analysis of data 
from multifragmentation, fission and compound nuclear decay experiments |48H51) . The values of these 
quantities obtained from n31o414 and n31o450 in Ref. [12] are displayed in Table ^ 






Figure 1: (Color online) Results for the free energy per nucleon F(T,p,S = 0), the pressure P{p,T,S = 0), the 
entropy per nucleon S{T, p,S = 0) and the internal energy per nucleon E{p, T,S = 0) in isospin-symmetric nuclear 
matter. The uncertainty bands correspond to calculations using two different sets of chiral low-momentum two- 
and three-body interactions, n31o414 (solid lines) and n31o450 (dash-dot lines). The unstable spinodal region is 
marked out explicitly. The critical point is shown as a circle (full circle for n31o414, open circle for n31o450). The 
zero-temperature endpoint of the low-density part of the spinodal is located at p ~ 2 • 10“^ fm“®. 


From the free energy per nucleon the pressure and the entropy per nucleon follow via standard ther¬ 
modynamic relations: 


P{T,p,S = 0)=p^ 


dF{T,p,6 

dp 


0 ) 


S{T,p,6 = 0) 


dF{T,p,6 = 0) 
dT 


( 5 ) 


The internal energy per nucleon is given hy E = F TS. The results for these quantities are shown in 
Fig. for temperatures in the range T = 0 — 25 MeV. The spinodal region^ where the homogeneous (i.e., 


^ Note that the value of the so-called critical compressibility factor is Zc = Pcl{Tc pc) — 0.29 for both n31o414 and n31o450; 
this is very similar to the values of Zc of various atomic or molecular fluids m, but differs from the value Zc = 0.375 

corresponding to equations of state of the van der Waals-Berthelot type m- _ 

^ In SNM the unstable spinodal region corresponds to {dP/dp)T ^ 0, with {dP/dp)T = 0 on the spinodal, cf. Sec. 
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single-phase constrained) system is unstable with respect to infinitesimal density fluctuations is shown 
explicitly. Note that for low temperatures the region of negative pressure extends into the metastable 
region (cf. also Fig. [I4| , which is a generic property of liquids that are self-bound at low temperatures 
and a well-known feature of superheated molecular liquids [54] . 



(MeV) 

Psat (fm“®) 

K (MeV) 

Tc (MeV) 

Pc (fm“®) 

Pc (MeV fm"®) 

n31o414 

-15.79 

0.171 

223 

17.4 

0.066 

0.33 

n31o450 

-15.50 

0.161 

244 

17.2 

0.064 

0.32 


Table I: Zero-temperature saturation energy Esat, density psat and incompressibility K, as well as 
the critical temperature Tc, density pc and pressure Pc of the liquid-gas phase transition in isospin- 
symmetric nuclear matter from the sets of chiral two- and three-body nuclear interactions n31o414 
and n31o450. 


III. PURE NEUTRON MATTER 

The free energy per particle in PNM, F{T,p,5 = 1), is obtained by restricting the isospin sum(s) in 
Eq. Q and in the different contributions in Eq. ([^. Moreover, in PNM the three-body contributions 
proportional to the low-energy constants ce,d and C4 are absent [38j . The results for the free energy per 
particle, the pressure, the entropy per particle and the internal energy per particle in PNM are shown 
in Pig. Note that the uncertainty bars obtained by varying the resolution scale (n31o414 vs. n31o450) 
increase with temperature and are significantly reduced as compared to the SNM results in Ref. [12j . 
which is due to the decreased magnitude of nuclear interactions in PNM. At very low temperatures the 
internal energy per particle increases monotonically with increasing density (as required by the absence 
of a liquid-gas instability in PNM), but otherwise there is a local minimum at finite density. 




P 




Figure 2: (Color online) Results for the free energy per particle F{T,p,S = 1), the pressure P{p,T,5 = 1), the 
entropy per particle S{T,p,S = 1) and the internal energy per particle E{p,T,S = 1) in pure neutron matter. 
The solid lines show the results from n31o414, the dash-dot lines the n31o450 results. The thick dashed lines 
correspond to the model-independent virial equation of state (VEoS) determined from neutron-neutron scattering 
phase shifts. The VEoS lines end where the fugacity is z = 0.5. 
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At very low energies, where higher partial waves are unimportant, the interaction between neutrons 
is characterized by the neutron-neutron scattering length as- In the regime where Os ~ — 19fm is large 
compared to the interparticle separation, 1 ^ [kpasl (with kp the Fermi momentum), a perturbative 
approach to neutron matter is not reliable. The model-independent virial equation of state (VEoS) 
computed by Horowitz and Schwenk in Ref. |55j from neutron-neutron scattering phase-shifts provides a 
benchmark for perturbative calculations of low-density neutron matter at nonzero temperature. In the 
virial expansion, the grand canonical expressions for the pressure and the density are expanded in powers 
of the fugacity 2 ; = exp(^/T), leading to 

P(r, z, <5 = 1) = ^ (z + z%iT) + Oiz^)) , p(T, z,6=l) = ^(^z + 2zH2(T) + 0{z^)) , (6) 

where /i is the chemical potential, and A = ^yj2TT)J{MT) is the neutron thermal wavelength. The second 
virial coefficient is given by 

^ 2 (r) = ^ 17 ^ dE eM-E/m StotiE) - 2-5/2, ( 7 ) 

where StotiE) is the sum of the isospin-triplet elastic scattering phase shifts at laboratory energy E. From 
the pressure and density as functions of the fugacity the free energy per particle F, entropy per particle 
S, and internal energy per particle E follow again from standard thermodynamic relations (see Ref. [55] 
for details). The results for these quantities are shown as green dashed lines in Fig. One sees that 
in the case of E, P and S there are almost no visible deviations between the VEoS and the perturbative 
results. This seemingly perfect agreement is however misleading, because the discrepancies corresponding 
to the different treatment of the interactions in the virial and the perturbative approach are overpowered 
by the large size of the (nonrelativistic) free Fermi gas contribution. The deviations are more transparent 
in the case of the internal energy per particle due to cancellations of the free Fermi gas terms in the free 
energy and entropy. The virial and perturbative results are closer at larger temperatures, since the EoS 
is less sensitive to the physics of large scattering lengths at higher momentum scales. 



P [fm 


Figure 3: (Color online) Interaction contribution to the internal energy per particle, Eint{p,T, 5 — 1), in pure 
neutron matter at T = lOMeV and low densities. The different lines correspond to the results from microscopic 
chiral nuclear interactions at first order (labeled “HF”) and second order in Eq. as well as the virial expansion 
truncated at second order (VEoS) and with uncertainty bands obtained by estimating the third-order term. The 
n31o414 and n31o450 results are almost identical. 

The differences between the virial and the perturbative results for the internal energy per particle 
are examined more closely in Fig. for T = lOMeV. To depict the deviations more clearly we have 


® Note that we have added the relativistic correction term to the VEoS lines. In particular, the zero-density limit of 

the internal energy per particle is given by E{T,p,5) Eq{T, p,5) |p—»o + p,&) |p—»o = 3T/2 -|- 15T2/(8M). 

This agrees with the expansion in powers of T of the internal energy per particle of a relativistic classical ideal gas, 
^classical = 3T + MKi{M/T)/K 2 {M/T) [where it’i, 2 (Vf/T) are modified Bessel functions]. 
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subtracted the noninteracting contributions, i.e., the quantity shown is i^int = E — Eq — i^rei- The 
virial results include uncertainty bands obtained from estimating the neglected third virial coefficient 
as |& 3 (T)| < \b2{T)\/2. We also show the perturbative results at the Hartree-Fock level [first order in 
Eq. ^]. One sees that compared to the Hartree-Fock results the inclusion of second-order contributions 
leads to much closer agreement with the virial expansion. The second-order calculation still slightly 
underpredicts the attractive interaction contributions, in contrast to the pseudopotential approach based 
on nucleon-nucleon scattering phase shift data that was explored in Ref. [S5]. We conclude that while 
the perturbative approach cannot fully capture the large scattering length physics of low-density neutron 
matter, the resulting errors are reasonably small when second-order contributions are included. 

In recent years, the zero-temperature EoS of PNM from chiral nuclear interactions has been studied 
by numerous authors within various many-body frameworks 01391 ESI ISlHMl- We compare our results 
to results obtained from perturbative calculations with various chiral interactions by the Darmstadt 
group (red band in Fig. 8 in Ref. [ST]) in Fig. In addition to the N2LO chiral three-neutron forces, 
their calculations also include all N3LO three- and four-neutron interactions. The uncertainty bands 
in their results were obtained by allowing large variations of the low-energy constants parameterizing 
the many-neutron forces. One sees that the (almost overlapping) results from n31o414 and n31o450 lie 
within these bands. In Fig. [^ we also show results obtained from auxiliary-field quantum Monte Carlo 
simulations with chiral N3LO two-body (AFQMC [NN]) and N3LO two-body plus N2LO three-body 
forces (AFQMC [NN-I-3N]) by Wlazlowski et al. [S^- The perturbative and the AFQMC results are very 
similar at densities p ^ 0.006fm~^, where both are in close agreement with the (fixed-node) quantum 
Monte Carlo calculations (based on the AV18 potential) of Gezerlis and Carlson [57]. However, at higher 
densities the EoS predicted by the AFQMC calculations (with three-body forces included) is significantly 
more repulsive. This discrepancy may be (partly) related to systematic errors in the AFQMC treatment 
(cf. also Ref. [55]). 



Figure 4: (Color online) Energy per particle in pure neutron matter at zero temperature, E{T = 0,p,S = 1), 
obtained from various many-body methods (see text for details). The inset magnifies the behavior at very low 
densities where quantum Monte Carlo simulations, labeled “QMC [AV18]”, are expected to be most accurate. 


IV. SYMMETRY FREE ENERGY, ENTROPY AND INTERNAL ENERGY 

From the results for the free energy per particle in homogeneous SNM and PNM the symmetry free 
energy Fsy^{T,p) is obtained via Eq. |H. The symmetry entropy and internal energy are related to the 
symnietry free energy via Ssym = —dEsym/dT and Egym = ^sym + TSgym. The results for Fgy^ TSsym 
and Egym are shown as functions of density at different temperatures in the left column of Fig. 0 In the 
insets we show the noninteracting contribution to these quantities, i.e., 

Fnonint,sym(T, p) = Fo{T, p, 1) - Fo{T, p, 0) -|- Frel{T, p, I) - F^el{T, p, 0), (8) 

in the case of the symmetry free energy. In the right column of Fig. [^we show Fsym(r,p), TSgy^{T,p), 
and Egy^{T,p) as functions of temperature at different densities. 
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Figure 5: (Color online) Left column: results for the symmetry free energy Fsyni{T, p), the symmetry entropy 
times temperature TSaymiT, p), and the symmetry internal energy Eaym{T,p), plotted as functions of density. 
The insets show the noninteracting (free Fermi gas) contribution to the different symmetry quantities. Right 
column: Fayin{T, p), TSaymiT, p), and Faym{T, p) as functions of temperature at different densities. The lines are 
interpolated, with calculated data points at T /MeV = 0, 3, 5,8,10,12,15, 20, 25. 

































































One sees that in the considered range of densities and temperatures, Fsym is a monotonic increasing 
function of density and temperature. The density and temperature dependence of T^sym is more involved. 
At low densities TSsym decreases monotonically with T, but for densities p > 0.2 fm~^ a local minimum 
is found at T ^ lOMeV. ^ The results from n31o414 and n31o450 are very similar for densities well 
below nuclear saturation density, but at higher densities the dependence on the resolution scale becomes 
significant. In particular, the decrease in the slope of Fsym with increasing density is more pronounced in 
the n31o450 results. The T dependence of TSsym{T, p) approximately balances that of F^ym{T,p), and as 
a result their sum, the symmetry internal energy E^ym, increases with density but varies only very little 
with temperature. At densities near nuclear saturation density the deviations of Esym(T,p ~ Psat) from 
its value at zero temperature are below 0.5 MeV. 

In Fig. we show the symmetry quantities with the noninteracting contributions subtracted, i.e., 

.Tint,sym(T, p) 

— ^sym {T,p) Tnonint,syni (T, p) , (9) 

as functions of temperature at different densities. In both cases the interaction contributions tend to 
counteract the temperature dependence of noninteracting contributions, cf. the insets in Fig. In the 
case of Fsym (and also TSsym) the noninteracing contributions dominate, but in the case of Esym the 
size of the noninteracting contributions and the ones from chiral nuclear interactions is more balanced, 
and the T dependence of both contributions approximately cancels each other, leading to the observed 
approximate temperature independence at densities near nuclear saturation density. 




Figure 6: (Color online) Temperature dependence of the interaction contributions to the symmetry free energy, 
Fsym,int(T, p), and the symmetry internal energy, Esyui,int{T, p) at different densities. The lines are interpolated, 
with calculated data points at T /MeV = 0, 3, 5,8,10,12,15, 20, 25. 


In Fig. we compare our results for the symmetry (free) energy at zero temperature to the results 
obtained by Drischler et al. m from calculations of the EoS of neutron-rich matter using several 
renormalization group-evolved chiral nuclear interactions. For comparison we also show the results from 
microscopic calculations within a variational approach by Akmal et al. based on the AVIS two-body 
and the Urbana UIX three-body potential. ^ While the results of Drischler et al. are compatible with our 
results, the calculations by Akmal et al. predict a symmetry energy that deviates visibly from the n31o414 
and n31o450 results. In Fig. we also show recent empirical constraints obtained from the analysis of 
isobaric analog states and neutron skins (lAS-l-NS) [69]. One sees that the n31o414 and n31o450 results 
lie in the lAS-l-NS bands in the entire constrained density region 0.04 < p/fm“^ < 0.16. 


We note that the temperature dependence of Fsym and TSsym approaches linear behavior in the limit of vanishing density, 
Fsym(T p—^ 0) = —TSsym(T P —^ 0) = F In 2. 

® The results by Akmal et al. include relativistic boost corrections as well as an artificial correction term added to reproduce 
the empirical saturation point of SNM (“corrected” in Table VI. and “A18-|-(5u-|-UIX* ” in Table VII. in Ref. [68]). 
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Figure 7: (Color online) Symmetry (free) energy as a function of density at zero temperature, Faym (T = 0,p). 
The results from the chiral nuclear interactions n31o414 and n31o450 are compared to those of Drischler et al. 
M and Akmal et al. [^. Also shown are empirical constraints from the analysis of isobaric analog states and 
neutron skins (lAS+NS). 


For densities close to nuclear saturation density the symmetry (free) energy at zero temperature is 
usually expanded around J = Fsym(T = 0,psat) in terms of a: = (p/psat — l)/3: 

Fsyni{T = 0, p) = J + Lx + —KsymX^ + 0{x^), (10) 

where L = SpsatdFsymiT = 0, p)/dp\p=p^^^ is called the slope parameter, and iCgym = 9pLt^^-^sym(r = 
0,p)/9p^ lp=PBat the symmetry incompressibility. The density where the ground state energy per particle 
in isospin-asymmetric nuclear matter has a local minimum is related to the parameters in the above 
expansion via Psat(^) — Psat[l — 3LS‘^/K] (cf. Ref. [70]).® The corresponding incompressibility K{S) 
obeys the approximate relation 

KiS)^K + KrS^, Kr=K,y^-6L, ( 11 ) 


where K^- is usually called the isobaric incompressiblity. In recent years, much effort has been invested 
in determining the parameters in Eqs. (101. The empirical values of J = 29.0 — 32.7MeV and to 
a lesser degree also L = 40.5 — 61.9MeV are relatively well constrained (values from [TT], see also 
lilMlIZHZl]), whereas experimental determinations of suffer from large uncertainties. For instance, 
from measurements of neutron skin thicknesses [T^ the value K^- = —SOOI^qq MeV was obtained, which 
is compatible with the giant monopole resonance measured in Sn isotopes m giving Kr = —550 ± 
100 MeV. Theoretical studies using a selection of Skyrme interactions however led to an estimate of 
Kt = —370 ± 120MeV [TD].^ Our results for J, L and ATt are given in Table [n| They are in agreement 
with the mentioned constraints. 



J (MeV) 

L (MeV) 

K.r (MeV) 

n31o414 

32.51 

53.8 

-424 

n31o450 

31.20 

48.2 

-434 


Table 11: Value of the (zero-temperature) symmetry energy at saturation density J, the slope 
parameter L, and the isobaric incompressibility Kt, extracted from the results obtained from the 
sets of chiral nuclear two- and three-body interactions n31o414 and n31o450. 


® The densities psat(<5) correspond to stable self-bound states only for isospin asymmetries up to the neutron drip point, 
<5 < <5nDi cf. Secs. [V^and [V^ 

^ However, in each case a slightly different definition of Kt is used, with the differences corresponding to higher-order terms 
in Eq. 0 and finite-size effects Oj. 
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V. THERMODYNAMICS OF ISOSPIN-ASYMMETRIC NUCLEAR MATTER 

In this section we examine the thermodynamic equation of state of isospin-asymmetric nuclear matter 
(ANM). The free energy per particle in homogeneous ANM is calculated as follows. The dependence of 
the nonrelativistic free Fermi gas contributions on the isospin asymmetry 6 is treated exactly, while the 
relativistic correction term and the interaction contributions are assumed to have a quadratic dependence 
on the isospin asymmetry: 


F(T, p, 6) ~ Fo(T, P, 6) + F,el(r, P, 0) + Fsym.rel(T, p) p, 0) + Fsym.int(T, p) 6^. (12) 

This approach is motivated in Sec. |V A[ where we investigate in detail the dependence on isospin asym¬ 
metry of the noninteracting contributions Fq{T,p, S) and F^g\(T, p, S). 

In Sec. |VB| we then discuss the construction of the spinodal in ANM and present our results for the 
trajectory of the critical temperature. The dependence of the neutron and proton chemical potentials 
on isospin asymmetry is examined in Sec. |V C[ and we determine the neutron drip point in cold nuclear 
matter. Finally, in Sec. EE we show results for the free energy per particle and pressure in isospin- 
asymmetric nuclear matter and determine the values of T and 6 where an isolated drop of liquid nuclear 
matter becomes unstable. 


A. Isospin dependence of free nucleon gas 

Here we examine the S dependence of the noninteracting contributions to the free energy per par¬ 
ticle in homogeneous ANM.® We compute the four leading terms in an expansion of Fq{T,p,6) and 
Fre\{T, p,S) in powers of 6^. The results show that the accuracy of the quadratic approximation for 
Fo{T,p,S) decreases significantly with increasing temperature, which necessitates the exact calculation 
of this contribution. The relativistic correction term on the other hand can be safely approximated via 
F^el{T, p, S) ~ Frel(T, P, 0) -f Fsym,rel{F^ p^ 6 . 

The noninteracting contributions to the free energy density, Fq = Fq + Fq and Fj-ei = F^l^ + can 
be expressed in terms of polylogarithms Li^(a;) = i.e., 

Fq^^{T, ln(-x„/p) Li3/2(a:n/p) - Li 5 / 2 (x„/p)) , (13) 

( 14 ) 

where are the neutron and proton effective one-body chemical potentials, Xn/p = — /T) 

and a = 2“^/^(M/7r)®/^. For given values of T, p and 5 the effective one-body chemical potentials 
are uniquely determined by p„/p = —ar®/^Li 3 / 2 (a;n/p)- 

The noninteracting contribution to the free energy per particle, Fjionint = {Fq F^f.\)/{Pn + Pp)) as a 
function of temperature T, nucleon density p and isospin asymmetry 5 can be expanded® in powers of 
around its value in SNM (<5 = 0): 

OO 

-^nonint (T, p, 5) = F„onint(T, p, 0) -f ^ F„onint.2n(T, p) 5^^. (15) 

n—1 


^ As noted in Sec. ^ in initial calculations we have found that at the densities and temperatures relevant for the liquid-gas 
phase transition the interaction contributions give rise to comparatively weaker terms with quartic and higher powers 
in 5. The contributions from two-nucleon interactions Fi nn F 2 nn (which give the dominant contribution to ^int 
at subnuclear densities) were found to be quadratic in 6 to high accuracy. In the case of -Fi, 3 n we have found that 
higher-order terms in <5 are more sizeable, but still small compared to those that emerge from ^0 finite T). Future 
research will quantify the isospin-asymmetry dependence of the interaction contributions in more detail. 

^ Odd-order terms in S in the expansion of -Fnonint arise only from the neutron-proton mass difference AM ~ 1.4 • 
which we neglect in this work. 
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The various expansion coefficients i?nonint, 2 n are given by 


B, 


nonint,2n 


iT,p) 


1 d^^FnonintiT,p,S) 

(2n)! 


(16) 


Setting ^ = 1 in Eq. 
coefficients, i.e., 


(15) we obtain the noninteracting symmetry free energy as the sum of the above 


-^nonint,sym(^;/^) ^ ^ 




(17) 


Here, we have introduced the weight factors 3nonint,2n = 7?nonint,2n/^nonint,sym as a means to specify the 
relative size of the different expansion coefficients. 

The rules of multivariable calculus lead to the following expression for the S derivative of order n of 
i G {0,rel}, at fixed density and temperature: 


^n^n/p\ 

dS^ ] 

^T,p 


(±l)”ar5/2(n- 1)! 

(1 ± d)" 


^^”^(a:„/p). 


(18) 


where the functions are defined recursively as 


^^"^(a:n/p) = 


f-n/p Li3/2(2:n/p) 9 


max(n-l,l) Lii/ 2 (a:„/p) dx^/^ 


'(x„/p)-(l-6„p)^^(”-^)(x„/p). 


n > 1, 

(19) 


with 6fc ; the Kronecker delta. The expressions to start the recursion are 

^ (s^n/p) = la(— Xjj/p) Li3/2(2:n/p) ~ Li5/2(a^n/p)i ^el (^n/p) = ~ ^^^Li7/2 (^^n/p) ■ (20) 

One then obtains for i?nonint. 2 n the expression 


Fnonint, 2 n(F, x) 




2n Li 3 / 2 (x) 


(2n) 


( 21 ) 


where x = — exp(/ro/r), with /tq the nucleon effective one-body chemical potential, which is uniquely 
determined by p = — 2 aT^/'^hi^/ 2 {x). 

The results for the first weight factor ( 3 nonint ,2 = H„onint.2/^nonint,sym and the ratios 
Pnonint.2n/Pnonint,2(n-i-i) for n = 1, 2, 3 Ere displayed in Fig. Note that the limits p —0 and T —>■ 0 do 
not commute. The p —^ 0 limit of the symmetry coefficients at finite temperature is given by 


Bnonint,2n{T 7^ 0, p —)■ 0) — 


T 5^" 
2(2n)!5(P 


((1 + (5) ln(l + 5) + (1 - d) ln(l - d)) 


T 


5=0 2n{2n — 1) ’ 


( 22 ) 


which comes entirely from the logarithmic terms, ~ ln(—x„/p) = in the expression for Fq = 

{Fq + FP)/(pn + Pp). The asymptotic behavior of the weight factors |3nonint,2r!, for T —>• oo is determined 
by the fact that the p —)■ 0 and T —>■ oo limits of |3nonint,2Tt = Hnonint, 2 n/Fnonint,sym coincide.As apparent 
from Fig. the asymptotic behavior sets in at relatively low values of T, causing the convergence rate of 
the quadratic expansion of ^)ionint to decrease significantly with increasing temperature. This behavior 
is entirely caused by the presence of the logarithmic terms in the nonrelativistic contribution, i.e., by 
the first term proportional to the sum of the effective one-body chemical potentials, ^ pp -I- pg. The 
p —>■ 0 limit of the second term in Fq proportional to [Li 5 / 2 (xn) -b Li 5 / 2 (xp)]/[Li 3 / 2 (xn) -b Li 3 / 2 (xp)] 
equals —T, which is independent of <5; the convergence rate of the expansion in powers of 6^ of this 
term alone increases very strongly with increasing temperature. Similarly, the p —>■ 0 limit of Fi-ei equals 


n/p 


This follows from the fact that x, 


—> 0 in both the p 0 and the T —> oo limit. 
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— 15T^/(8M), and the quadratic approximation of the relativistic correction term becomes increasingly 
accurate with increasing temperature. 

For comparison, in Fig.[^we show the results for the weight factors anonint. 2 n in the analogous expansion 
in powers of of the noninteracting contributions to the internal energy per particle, Eq and Erei- The 
noninteracting contributions to the neutron and proton internal energy densities are given by 




3aT5/2 


^^5/2 (^n/p) 7 


75aT7/2 

16M 


Li7/2(2;n/p) ' 


Ti5/2(^n/p) Ti3/2(^n/p) 
16M Lii/2(a;n/p) 


(23) 

(24) 


Both Eq —)■ 3T/2 and E^ei 15T^/(8M) become independent of (5 as p —>■ 0, and the convergence rate 
of the expansion of both terms increases strongly with increasing temperature (the increase is signifi¬ 
cantly more pronounced in the case of Eq, which is reflected in the increase of the deviations between 
relativistically improved and the nonrelativistic results with increasing values of n). 



T [MeV] 





0 5 10 15 20 25 

T [MeV] 


Figure 8: (Color online) Temperature dependence of the first weight factor |3nonint,2 and the ratios 
Pnonint, 2 n/Pnonint, 2 (n+i) for n = 1,2,3 at different densities, corresponding to the expansion of Fnonint{T, p, 5) 
in powers of S^. The full lines show the results with the relativistic correction term included, the dotted lines 
the nonrelativistic results. Note that the deviations between the relativistically improved and the nonrelativistic 
results decrease with increasing values of n, indicating the opposite convergence behavior of the expansion in 
powers of 5^ of Fq and Flei, respectively. 
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Figure 9: (Color online) Temperature dependence of the first weight factor anonint ,2 and the ratio 
«nonint,2/anonint,4, Corresponding to the expansion in powers of 5^ of the noninteracting contributions to the 
internal energy per particle, Eq and Srei- Note the logarithmic scale in the second plot. 


B. Spinodal and critical temperature 


The onset of spinodal instability is associated with the violation of a number of equivalent stability 
criteria [771 [TS], which are derived from the fundamental principle of maximum entropy. In the canonical 
representation the corresponding stability requirement is that the free energy density F = pF at fixed 
temperature T is a convex function of the component densities pi and p 2 ^ implying that the Hessian 
matrix Fij has no negative eigenvalues. In the unstable region delineated by the spinodal the (analytical) 
free energy density is concave; in the metastable region between the spinodal and binodal (coexistence 
boundary) the free energy density is locally convex and the system is protected against phase separation 
by a nucleation barrier. The Hessian matrix Fij is given by 


pi, P 2 ) 


ld^F{Fp,,p2)] 


'dpt{T, Pi,P2y 

dpidpj 


dPj 


biG{l,2}. 


(25) 


Its eigenvalues are given by 


i±{T,pi,p2) 

tr[J-y] ± (tr[J-,,f - 4Aet[F^j]f'^ 


1 

“2 

Fu + F 22 ± ((-Fn - F 22 f + 47^12 

^1/2' 


(26) 


The signs of the eigenvalues are invariant under (linear) basis transformations. From the data F{T,p,S) 
they are readily evaluated using as independent density parameters pi = Pn + Pp — P (nucleon density) 
and P 2 = Pn~Pp = pS (isospin asymmetry density). In this basis the Hessian matrix becomes diagonal at 
6 = 0 with eigenvalues = {d'^F/dpl)T,pi\p 2 =o > 0 and = (i9^-F/i9p^)T.p2 Ip 2 =o = P~^idP/dp)T,s\6=o 
(this result depends on the neglect of isospin-symmetry breaking effects). Hence, in SNM the region 
inside the spinodal corresponds to a negative isothermal compressibility kt = p~^{dp/dP)T,s, where 
> 0 is a stability criterion for a pure substance. 

The exact expressions for the nonrelativistic free Fermi gas contribution to the Hessian matrix compo¬ 
nents are given by^^ 


-^0,11(2^) Moi Po) 
Po. 12{T, PoyPo) 


T-1/2 / I 

4a VLii/ 2 (a:n) 

T-1/2 / I 

4a VLii/ 2 (a:n) 


1 

Lii/2(x 

1 

Lii/2(a:F 


— Po,22(P, PoFo)j 


(27) 

(28) 


where again x^/p = — exp(p,g'^^/T) and a = 2 ^/^(M/tt)^/^. In the limit i5 —1 the proton density van¬ 
ishes, Pp 0, and the proton effective one-body chemical potential diverges (at finite T), /Tq —>■ — 00 , thus 


11 


Note that Jb.ll = Jb,22, but Tint, 11 i=- -^int,22- 
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Lii/ 2 (a^p) —>■ 0- Hence, the exact calculation of the nonrelativistic free Fermi gas contribution Fq{T,p,S) 
leads to divergent behavior of the Hessian components in the limit of vanishing proton concentration. The 
same divergent behavior is obtained in the exact calculation of at zero temperature. The unstable 
region then vanishes at a value 5 < 1 for all values of T. This constraint is lost if the free Fermi gas 
contribution is approximated by truncating the expansion in powers of <5^ of Fq(T,p,S) at a finite order, 
e.g., at first order as in the usual quadratic isospin approximation. 

The evolution of the critical temperature Tc{d) where (for a given value of (5) the unstable concave 
region vanishes is depicted in Fig. The results from n31o414 and n31o450 are very similar; in both 
cases the critical lines end approximately at an isospin asymmetry 5®”'^ ~ 0.9994 or a proton concentration 


yp®"®^ ~ 3 • 10“^. The value 5®“'^ ~ 0.9994 exceeds the critical line endpoints obtained in Refs. [TSlFMfSO] 
using different phenomenological models. We note that in nuclear matter with 6^0 the coexistence 
region does not vanish at the critical temperature Tc{S) but at a higher temperature I^x(<5), the so- 
called maximum temperature umig. The existence of a neutron drip point (see Sec. |V C[ ) entails that 
at zero temperature the binodal extends to <5 = 1 over a finite region of densities or pressures. The 
trajectory of the maximum temperatures therefore reaches its zero-temperature endpoint at vanishing 
proton fraction = 0. 

For comparison, in Fig. 10 we also show the trajectories of the temperature {S) where the region with 


negative isothermal compressibility vanishes at fixed S. The exact expression for the nonrelativistic 


free Fermi gas contribution to is given by 




^T,0 ~ 


(Li3/2(a^n) + Li3/2(a;p)) 


(1 + ^)^ 

Lii/2(a;n) 


Lii/2(a;p) 


(29) 


For both n31o414 and n31o450 the Tkt(^) trajectories end at approximately <5®"'^ ~ 0.82 or — 0.09, 
which exceeds the values Fp® — 0-053 from Ref. [ST] obtained in an in-medium chiral perturbation 
approach and Fp — 0.045 from Ref. [5S] obtained by applying the functional renormalization group to 
a chiral nucleon-meson model. 



Figure 10: Trajectories of the critical temperature Tc{S) determined from n31o450 and n31o414. The trajectories 
end at (5 ~ 0.9994. Also shown are the trajectories of the temperatures TKr^{S) where the region with negative 
isothermal compressibility kt vanishes. The calculated data points are shown explicitly. 


C. Stable self-bound liquid 


From the data F(T,p,S) the neutron and proton chemical potentials are obtained via 

,, dF{T,p,S) , It SdF{T,p,S) 

fJ-n/p[T,p,d) = --± 


dp 


dS 


(30) 


The results for Pn{T, p, 5) and Pp{T, p, 6) are displayed in Fig. f^for temperatures T /MeV = 0,15. One 
sees that Pn{T, p, 6) increases and Pp{T, p, i5) decreases with <5. Tlie chemical potentials at finite T diverge 
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as p —>■ 0, but at zero temperature p„/p —>■ 0 for p —0. The origin of this feature is the collapse of Fermi- 
Dirac distribution functions into Heaviside step functions at T = 0, which eliminates the logarithmic 
divergence of the free energy per particle at vanishing density [HIIS]. 

The Gibbs conditions for the coexistence of two bulk phases I (liquid) and II (gas) in mutual thermo¬ 
dynamic equilibrium are 


yl ^ rj.ll 


A, = m", Mp = fj^p- 


(31) 


Whereas at finite T liquid-gas equilibrium corresponds to finite values of density p and proton fraction 
Y in both phases, the vanishing of Pn/p(T = 0,p,S) at vanishing density entails that at T = 0 the Gibbs 
conditions for the neutron and proton chemical potentials can in most cases not be satisfied, leading to a 
gas phase that is either empty (vacuum) or contains only neutrons [liiiHa. The neutron drip point, ^nd, 
is given by the value of S where the neutron chemical potential at vanishing temperature and pressure 
becomes positive. For isospin asymmetries S < Jnd an isolated drop of cold liquid nuclear matter is 
stable (in equilibrium with the vacuum), defining a stable self-bound state. As seen from Fig. 11 in our 


equation of state neutron drip occurs at an isospin asymmetry 5 nd — 0.30 or a proton concentration 
Tp_ND = (1 ~ <5nd)/2 — 0.35, which is similar to results obtained with effective Skyrme interactions 

MM. 



Figure 11: (Color online) Neutron and proton chemical potentials, pn{T,p,5) and fip{T, p,S), in (homogeneous) 
isospin-asymmetric nuclear matter at temperatures T/MeV = 0,15, calculated using n31o414 (solid lines) and 
n31o450 (dash-dot lines). Stable self-bound states are shown as thick bright red lines with circles (full circles for 
n31o414, open circles for n31o450); the lines end at the neutron drip point 5 nd — 0.30. 
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D. Metastable self-bound liquid 


The zero-temperature results for the (ground state) energy per particle E = F and the pressure 
P = p^dF/dp are displayed in Fig. 12 as functions of the nucleon density p for different values of 5. 


The trajectory of the points where the energy per particle has a local minimum and the pressure is zero 
is shown explicitly. These points correspond to the properties of a drop of cold liquid nuclear matter 
surrounded by vacuum. For 5 < 5nd — 0.30 the cold drop is stable (the local energy minimum lies on 
the binodal, cf. Sec. VC I, and for Jnd < S < (5fp — 0.66 the local energy minimum lies in the metastable 
region between the binodal and the spinodal. We refer to the point <5fp — 0.66 where the trajectory of the 
local energy minima encounters the spinodal as the fragmentation point (FP). The energy per particle 
at neutron drip and at the fragmentation point is i^ND — —IS.OMeV and Epp ~ —3.1 MeV, respectively. 
For comparison we follow the local energy minima also into the unstable spinodal region; i.e., we show 
also the point where both derivatives of the (analytical) energy per particle vanish (saddle point, SP) 
at 5sp — 0-81 as well as the local energy minimum at (5 ~ 0.76. For S > 0.76 the energy per particle is 
positive at all (finite) densities, and for S > Jsp the pressure is a semipositive definite function of density. 



Figure 12: (Color online) Energy per particle E = F and pressure P in homogeneous isospin-asymmetric nuclear 
matter at zero temperature, calculated using n31o414 (solid lines) and n31o450 (dash-dot lines). The trajectories 
of the local energy minima are shown as thick dark gray (bright red below neutron drip) lines with circles (full 
circles for n31o414, open circles for n31o450). The trajectories end at the fragmentation point (5 fp — 0.66. The 
inset magnifies the behavior of the pressure at low densities. 


The finite-temperature results for the free energy per particle F{T, p, i5) and the pressure P{T, p, (5) are 
shown in Fig. 13 for temperatures T/MeV = 5,15. At finite T the trajectories of the local free energy 
minima lie entirely in the metastable region; an isolated drop of hot liquid nuclear matter has to be 
stabilized by a surrounding nucleon gas. At T = 5 MeV the trajectory ends at Jfp — 0.61, defining the 
fragmentation temperature for nuclear matter with proton concentration Yp ~ 0.195. The T = 5 MeV 
saddle point is located at Jsp — 0.68 for n31o414 and 5 sp — 0.69 for n31o450. For T = 15 MeV no local 
free energy minimum exists; for T > 13.5 MeV the analytical free energy per particle is a monotonic 
increasing function of density for all values of 6 (cf. Fig. 151. 

The relation between the spinodal, the binodal, and the trajectories of the local free energy minima 
is illustrated in Fig. The two plots in Fig. 14 represent isoplethal {S = const) and isothermal cross 
sections of the respective surfaces (spinodal, binodal, surface of local free energy minima) in (T, p, S) 
space (cf. also Refs. (TS] HZl)- In the second plot we also show the surface with divergent isothermal 
compressibility kt = p~^{dp/dP)T,s, which corresponds to the violation of the stability criterion > 0 
for a one-component system. Hence, the saddle point (SP) where both derivatives of the free energy per 
particle with respect to the nucleon density vanish coincides with the fragmentation point (FP) where 
a liquid drop becomes unstable only for (5 = 0 where nuclear matter behaves like a pure substance. For 
S ^ 0 the more restrictive two-component stability criteria are needed m ( > 0 is not a relevant 

stability criterion in that case) and the SP is located in the interior of the spinodal. 

Finally, in Fig. 15 the trajectory of the fragmentation temperatures Tfp(( 5) is shown [for comparison we 


also show the saddle point temperatures rsp(5)]. This trajectory determines the range of temperatures 
and isospin asymmetries for which a self-bound liquid state exists. 
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Figure 13: (Color online) Free energy per particle F{T, p, S) and pressure P{T, p, 5) at temperatures T = 5,15 MeV. 
The solid lines show the n31o414 results, the dash-dot lines the n31o450 results. The thick dark lines with circles 
depict the trajectories of the local free energy minima, up to the point where they encounter the spinodal. The 
insets magnify the behavior of the pressure at low densities. 




Figure 14: (Color online) Left plot: binodal, spinodal, and trajectory of local free energy minima (“self-bound 
liquid”) in SNM (5 = 0). Right plot: T — 10 MeV cross sections of the spinodal, the kt = oo boundary, and the 
surface of local free energy minima (“self-bound liquid”); only the 5 = 0 endpoints of the binodal are shown (we 
did not construct the binodal for 5 yf 0). The critical points (CP), fragmentation points (FP) and saddle points 
(SP) are shown explicitly in both plots (in SNM the FP and the SP coincide). 
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Figure 15: (Color online) Trajectory of the fragmentation temperature Tfp{S) above which no metastable self¬ 
bound state can exist (lower red line). For comparison we also show the trajectory of the points where the 
(analytical) free energy per particle has a saddle point, Tsp{5) (upper blue line). The calculated data points are 
shown explicitly. 


VI. SUMMARY 

In this work, we have investigated in detail the temperature and density dependence of the symmetry 
free energy Fsyni{T,p) = F{T,p,6 = 1) — F{T,p,6 = 0) in homogeneous nuclear matter using chiral 
effective field theory interactions constructed at resolution scales A = 414,450 MeV. The free energy per 
particle of isospin-symmetric nuclear matter, F{T,p,S = 0), and of pure neutron matter, F(T,p,S = 1), 
have been calculated in second-order many-body perturbation theory (Kohn-Luttinger-Ward formalism). 
Constraints from the nuclear saturation point, the critical point of the liquid-gas phase transition, and 
the density dependence of at zero temperature are reproduced, and our results are in reasonable 
agreement with the virial expansion of the neutron matter equation of state at low fugacities. 

The four leading coefiicients in an expansion of the noninteracting contributions Fq and J^rei in terms 
of 6^ have been examined, and we have found that the convergence rate of the expansion of decreases 
significantly with temperature. Therefore we have used the exact expressions for Fq in computing the 
free energy per particle F(T, p, S) in isospin-asymmetric nuclear matter. The many-body contributions 
from nuclear interactions on the other hand have been assumed to have a quadratic dependence on the 
isospin asymmetry S. 

From the results for F{T^ p, <5) we have computed the pressure P{T, p, <5) and the neutron and proton 
chemical potentials /in/p(T, p, d), and we have constructed the trajectory of the critical temperature Tc{S). 
The critical line ends at a proton fraction ~ 3 • 10“^. The neutron drip point in infinite nuclear 

matter at zero temperature has been located at Yp^ND — 0.35. Furthermore, we have determined the 
trajectory of the fragmentation temperature Tfp( 5) above which no metastable self-bound state exists. 

Future work will be aimed at improving the description of the isospin-asymmetry dependence of the 
interaction contributions to the thermodynamic equation of state, a more detailed treatment of the low- 
density region including Coulomb and surface effects, and the extrapolation of the equation of state to 
higher temperatures and densities required for simulations of core-collapse supernovae and binary neutron 
star mergers. 
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